---
title: "Civil-Military Relations in the Aftermath of Coups Replication Code"
author: "Artem Kyzym"
date: Updated on 24/10/2025
output:
  pdf_document: 
    latex_engine: xelatex
  word_document: default
  html_document:
    df_print: paged
---

LOAD DATA

```{r}
options(scipen = 999)

library(pacman)
p_load(tidyverse, foreign, haven, MASS, countrycode, plm, 
       miceadds, skimr, devtools, pdynmc, stats, miceRanger, 
       did, did2s, showtext, panelr, ggpubr, modelsummary, 
       pandoc, data.table, kableExtra, htmltools, styler)

colpus_cw <- read.csv("colpus_cw_final.csv")
pt_cw <- read.csv("pt_cw_final.csv")
```

MASTER THEME

```{r}
showtext_auto()
au_theme <- theme_bw() +
  theme(plot.title = element_text(color = "#24496B", family = "sans",
                                  hjust = 0.5, size = 41),
        axis.title.x = element_text(hjust = 0.5, size = 37,
                                  color = "#24496A", family = "sans"),
        axis.title.y = element_text(hjust = 0.5, size = 37,
                                    color = "#24496A", family = "sans"),
        axis.text.x = element_text(family = "sans", color = "#24496A"),
        axis.text.y = element_text(family = "sans", color = "#24496A"),
        panel.border = element_blank(),
        axis.ticks.y = element_blank(),
        axis.ticks.x = element_blank(),
        plot.caption = element_text(color = "#838383"),
        legend.position = "right",
        text = element_text(size = rel(7)),
        strip.text.x = element_text(size = rel(6)),
        strip.text.y = element_text(size = rel(6)),
        legend.text = element_text(size = rel(4.5), family = "sans"),
        legend.title = element_text(size = rel(5.25), family = "sans"))

stderr <- function(x, na.rm = FALSE) {
  if (na.rm) x <- na.omit(x)
  sqrt(var(x) / length(x))}
```

DATA PRE-PROCESSING

```{r}
#Random Forest Imputation for Covariates
colpus_vectors <- subset(colpus_cw, select = -c(1:14, 47:56))

miceObj <- miceRanger(colpus_vectors, m = 5, returnModels = T,
                      verbose = T, seed = 777)

colpus_vectors_new <- impute(colpus_vectors, miceObj, verbose = T)
colpus_vectors_final <- colpus_vectors_new$imputedData$Dataset_1

#Diagnostics
plotCorrelations(miceObj)
plotDistributions(miceObj)
plotImputationVariance(miceObj)
plotModelError(miceObj)
plotVarImportance(miceObj)

#Merge back With Original Variables
colpus_dep <- subset(colpus_cw, select = c(1:14, 47:56))
colpus_cwrf <- as.data.frame(cbind(colpus_dep, colpus_vectors_final))

#Factor Fixed Effects
colpus_cwrf$leaderspell <- as.factor(colpus_cwrf$leaderspell)
colpus_cwrf$year <- as.factor(colpus_cwrf$year)

#Counterbalancing (Sample Average = 1.468)
colpus_cwrf <- colpus_cwrf %>%
  mutate(cb_count_1 = cbcount + 0.01) %>%
  filter(year %in% c(1960:2010)) %>%
  mutate(mean_cb = mean(cb_count_1, na.rm = T)) %>%
  mutate(counterbalancing = ifelse(cb_count_1 > mean_cb, 1, 0))

#Log Counterweights
colpus_cwrf <- colpus_cwrf %>% mutate(cb_log = log(cbcount + 1))

#Past Coups at Regime Level
colpus_cwrf <- colpus_cwrf %>%
  group_by(leaderspell) %>%
  mutate(past_coups = ifelse(post_fail_colpus == 0, max(sum_any), 0))
```

DESCRIPTIVE STATISTICS (FIGURE 1 AND FIGURE 2)

```{r}
#Counterweights by Year
colpus_cwrf_desc <- colpus_cwrf %>%
  group_by(year) %>%
  filter(year %in% c(1960:2010)) %>%
  dplyr::summarise(mean_cb = mean(cbcount, na.rm = T))

ggplot(colpus_cwrf_desc, aes(x = year, y = mean_cb)) +
  geom_line(aes(group = 1), linewidth = 1) + au_theme +
  geom_ribbon(aes(group = 1, ymax = mean_cb, 
                  ymin = 1), alpha = 0.1) +
  labs(y = "Mean Number of Counterweights") +
  theme(axis.title.x = element_blank()) +
  coord_cartesian(ylim = c(1, 2)) +
  scale_y_continuous(breaks = seq(1, 2, 0.25)) +
  scale_x_discrete(breaks = seq(1960, 2010, 5)) +
  geom_hline(yintercept = 1.468, lty = 2)
#ggsave("timetrend2.png", width = 8, height = 4.5, dpi = 300)

#Distribution of Failed Coups
colpus_coups <- colpus_cwrf %>%
  filter(colpus_fail == 1) %>%
  group_by(year) %>%
  dplyr::summarise(count_coup = n())

colpus_none <- as.data.frame(c(1986, 1987, 1999, 2005, 2007, 2008, 2009))
colnames(colpus_none) <- c("year")
colpus_none$count_coup <- 0

colpus_joint <- as.data.frame(rbind(colpus_coups, colpus_none))
colpus_joint <- colpus_joint %>% arrange(year)

ggplot(colpus_joint, aes(x = year, y = count_coup)) +
    geom_segment(aes(x = year, xend = year, y = 0, 
                     yend = count_coup), color = "black") +
  geom_point(color = "black") + au_theme +
  labs(y = "Number of Failed Coups") +
  theme(axis.title.x = element_blank()) +
  coord_cartesian(ylim = c(0, 8)) +
  scale_y_continuous(breaks = seq(0, 8, 1)) +
  scale_x_discrete(breaks = seq(1950, 2010, 10))
#ggsave("countcoup.png", width = 8, height = 4.5, dpi = 300)

#Average Leader-Spell Length (9.05 Years)
colpus_leader <- colpus_cwrf %>%
  group_by(leaderspell) %>%
  summarise(number = n())
mean(colpus_leader$number)

colpus_desc <- colpus_cwrf %>%
  ungroup() %>%
  dplyr::select(c(
    counterbalancing, one_cw, cb_log, t_leader,
    latent_personalism, rgdppc, g_growth_rate,
    oil_new, per_solider, mid, acd_intrastate,
    exclpop, milethnic_hetero, defense))

colpus_desc <- skim(colpus_desc)
colpus_desc %>% filter(skim_type == "numeric")

sum(colpus_cwrf$colpus_fail, na.rm = T)
unique(colpus_cwrf$leaderspell) #Leader-spells
unique(colpus_cwrf$gwf_casename) #Regimes
```

DEFINE TREATMENT (FIGURE 3)

```{r}
#Never-Treated Control (N = 2879)
colpus_controls <- colpus_cwrf %>%
  filter(colpus_fail == 1) %>%
  pull(leaderspell)

colpus_control <- colpus_cwrf %>% 
  filter(!leaderspell %in% c(colpus_controls))
colpus_control$never <- "Never-Treated"

#Average Length = 8.59 Years
colpus_control_length <-
  colpus_control %>%
  group_by(leaderspell) %>%
  summarise(number = n()) 
mean(colpus_control_length$number) 

#Untreated (N = 296)
colpus_treatment <- colpus_cwrf %>% 
  filter(leaderspell %in% c(colpus_controls))

colpus_untreated <- colpus_treatment %>% 
  filter(post_fail_colpus == 0)
colpus_untreated$never <- "Untreated"

#Average Length = 5.02 Years
colpus_untreated_length <-
  colpus_untreated %>%
  group_by(leaderspell) %>%
  summarise(number = n()) 
mean(colpus_untreated_length$number)

#Treated (N = 724)
colpus_treated <- colpus_treatment %>% 
  filter(post_fail_colpus == 1)
colpus_treated$never <- "Treated"

#Average Length = 7.7 Years
colpus_treated_length <-
  colpus_treated %>%
  group_by(leaderspell) %>%
  summarise(number = n())
mean(colpus_treated_length$number) 

#Mean CW By Group
colpus_regroup <- as.data.frame(rbind(
  colpus_control, colpus_untreated, colpus_treated))

mean_group_weights <- 
  colpus_regroup %>%
  group_by(never) %>%
  dplyr::summarise(
    mean_cw = mean(cbcount, na.rm = T),
    se_cw = stderr(cbcount, na.rm = T))

mean_group_weights$low <- mean_group_weights$mean_cw - (1.96 * mean_group_weights$se_cw)
mean_group_weights$high <- mean_group_weights$mean_cw + (1.96 * mean_group_weights$se_cw)
mean_group_weights$order <- as.character(c(1, 3, 2))
mean_weights_labels <- c("Never-Treated", "Untreated", "Treated")

ggplot(mean_group_weights, aes(x = order)) +
  geom_pointrange(aes(y = mean_cw, ymin = low, ymax = high)) +
  au_theme + coord_cartesian(ylim = c(1, 2)) +
  scale_y_continuous(breaks = seq(1, 2, 0.2)) +
  scale_x_discrete(labels = mean_weights_labels) +
  labs(y = "Mean Number of Counterweights", x = "Treatment Status",
       caption = "Note: 95% Confidence Intervals Shown") +
  theme(plot.caption = element_text(size = 30, hjust = 0))
#ggsave("mean_group.png", width = 8, height = 4.5)

g1 <- subset(mean_group_weights, order == "3")

ggplot(mean_group_weights, aes(x = order)) +
  geom_pointrange(aes(y = mean_cw, ymin = low, ymax = high)) +
  geom_pointrange(data = g1, aes(y = mean_cw, ymin = low, 
                                 ymax = high), color = "#C30010") +
  au_theme + coord_cartesian(ylim = c(1, 2)) +
  scale_y_continuous(breaks = seq(1, 2, 0.2)) +
  scale_x_discrete(labels = mean_weights_labels) +
  labs(y = "Mean Number of Counterweights",
       caption = "Note: 95% Confidence Intervals Shown") +
  theme(plot.caption = element_text(size = 30, hjust = 0),
        axis.title.x = element_blank())
#ggsave("mean_group2.png", width = 8, height = 4.5, dpi = 300)

#T-test (T-value of 2.7 (p < 0.01))
sd(colpus_treated$cbcount, na.rm = T)
sd(colpus_untreated$cbcount, na.rm = T)
(1.778431 - 1.422111) / ((1.195902 / sqrt(724)) + (1.505114 / sqrt(296))) 
```

DESCRIPTIVE PARALLEL TRENDS (FIGURE 4)

```{r}
#Synthetic Control For All Failed Coups
slimdata <- colpus_cwrf %>% dplyr::select(
  gwf_country, ccode, leaderspell, year, t_leader,
  colpus_fail, post_fail_colpus, cbcount)

slimdata <- slimdata %>%
  group_by(leaderspell, t_leader) %>%
  mutate(ref_yr = case_when(
    colpus_fail == 1 ~ t_leader,
    TRUE ~ NA_integer_)) %>%
  ungroup() %>%
  group_by(leaderspell) %>%
  fill(ref_yr, .direction = "downup") %>%
  mutate(yr_diff = abs(ref_yr - t_leader))

slimdata_avg <- slimdata %>%
  filter(post_fail_colpus == 0) %>%
  mutate(yr_diff = yr_diff - (yr_diff * 2)) %>%
  dplyr::select(c(ccode, year, yr_diff))

slimdata <- left_join(slimdata, slimdata_avg,
  by = c("ccode", "year", "leaderspell"))

slimdata <- slimdata %>%
  mutate(duration = if_else(is.na(yr_diff.y), yr_diff.x, yr_diff.y)) %>%
  #filter(complete.cases(duration)) %>% filter(duration %in% c(-10:10))
  mutate(pool = if_else(is.na(duration), "Donor", "Treated"))

slimdata <- slimdata %>%
  group_by(duration) %>%
  mutate(
    before_after_treated = mean(cbcount, na.rm = T),
    before_after_sd = stderr(cbcount, na.rm = T))

placeboframe <- tibble(
  time = double(),
  counterbalancing = double(),
  sd_counterbalancing = double(),
  potentialcoupyear = double())

for (c_year in 1960:2010) {
  print(c_year)
  newdata <- slimdata %>%
    filter(pool == "Donor") %>%
    mutate(newyear = as.numeric(as.character(year)),
           time = newyear - c_year) %>%
    filter(time > -10 & time <= 10) %>%
    group_by(time) %>%
    dplyr::summarize(
      counterbalancing = mean(cbcount, na.rm = TRUE),
      sd_counterbalancing = sd(cbcount, na.rm = TRUE) / sqrt(n())) %>%
    mutate(potentialcoupyear = c_year)
  placeboframe <- rbind(placeboframe, newdata)
}

treated_data <- slimdata %>%
  filter(pool == "Treated") %>%
  group_by(leaderspell) %>%
  mutate(coup_year = ifelse(any(duration == 0),
    as.character(year[duration == 0]), NA)) %>%
  filter(duration > -10 & duration <= 10)

colnames(placeboframe) <- c("duration", "donor_cb", "donor_sd", "coup_year")
placeboframe$coup_year <- as.character(placeboframe$coup_year)

megaframe <- left_join(treated_data, placeboframe, 
                       by = c("duration", "coup_year"))

#Summarize by Everything
final_frame <- megaframe %>%
  group_by(duration) %>%
  dplyr::summarize(
    treated_mean = mean(before_after_treated, na.rm = T),
    treated_sd = mean(before_after_sd, na.rm = T),
    control_mean = mean(donor_cb, na.rm = T),
    control_sd = mean(donor_sd, na.rm = T))

treated_trends <- final_frame %>% 
  dplyr::select(c(duration, treated_mean))
treated_trends$status <- "1Treated"
colnames(treated_trends) <- c("duration", "mean_cw", "status")

control_trends <- final_frame %>%
  dplyr::select(c(duration, control_mean))
control_trends$status <- "2Never-Treated"
colnames(control_trends) <- c("duration", "mean_cw", "status")

parallel_trends <- as.data.frame(rbind(treated_trends, control_trends))
labss <- c("Treatment", "Never-Treated")

ggplot(parallel_trends, aes(x = duration, y = mean_cw, group = status)) +
  geom_line(aes(color = status), linewidth = 0.5) + au_theme +
  geom_point(aes(color = status, linewidth = 1)) +
  labs(y = "Mean Number of Counterweights",
       x = "Number of Years Before and After a Failed Coup",
       color = "Group") +
  scale_color_manual(values = c("#C30010", "darkblue"), label = labss) +
  scale_x_continuous(breaks = seq(-9, 10, 1)) +
  scale_y_continuous(breaks = seq(1.0, 2.0, 0.1)) +
  geom_vline(xintercept = 0, linetype = 2) +
  theme(legend.position = "bottom")
#ggsave("paralleltrend11.png", width = 8, height = 4.5, dpi = 300)

ggplot(parallel_trends, aes(x = duration, y = mean_cw, group = status)) +
  geom_smooth(method = "loess", se = T, level = 0.95, span = 0.5,
              aes(color = status), size = 1) + au_theme +
  labs(y = "Mean Number of Counterweights",
       x = "Number of Years Before and After a Failed Coup",
       color = "Group") +
  scale_color_manual(values = c("#C30010", "darkblue"), label = labss) +
  scale_x_continuous(breaks = seq(-9, 10, 1)) +
  scale_y_continuous(breaks = seq(1.0, 2.0, 0.1)) +
  geom_vline(xintercept = 0, linetype = 2) +
  theme(legend.position = "bottom")
#ggsave("paralleltrend2.png", width = 8, height = 4.5)
```

MAIN RESULTS (TABLE 1)

```{r}
#Time-To-Treatment
before_after_es <- colpus_cwrf %>%
  group_by(leadername, t_leader) %>%
  mutate(ref_yr = case_when(
    colpus_fail == 1 ~ t_leader,
    TRUE ~ NA_integer_)) %>%
  ungroup() %>%
  group_by(leaderspell) %>%
  fill(ref_yr, .direction = "downup") %>%
  mutate(yr_diff = abs(ref_yr - t_leader))

before <- before_after_es %>%
  filter(post_fail_colpus == 0) %>%
  mutate(yr_diff = yr_diff - (yr_diff * 2)) %>%
  dplyr::select(c(ccode, year, yr_diff))

before_after_es <- left_join(before_after_es, before,
  by = c("ccode", "year", "leaderspell"))

before_after_es <- before_after_es %>%
  mutate(duration = if_else(is.na(yr_diff.y), yr_diff.x, yr_diff.y))
before_after_es$duration[is.na(before_after_es$duration)] <- -Inf
before_after_es <- before_after_es %>% arrange(ccode, year)

#Anticipation
before_after_es <- before_after_es %>%
  mutate(anticipation = ifelse(duration == -1, 1, 0),
         anticipation2 = ifelse(duration == -2, 1, 0))

#Time Since Last Successful Coup
before_after_es <- before_after_es %>%
  arrange(ccode, year) %>%
  group_by(ccode) %>%
  mutate(
    time_since_last_coup = {
      counter <- 0
      sapply(1:n(), function(i) {
        if (colpus_success[i] == 1) {
          counter <<- 1
        } else {
          counter <<- counter + 1
        }
        counter
      })
    }
  ) %>%
  ungroup()

#Lag Covariates Within Leader-Spells
dt_colpus <- as.data.table(before_after_es)
dt_colpus <- dt_colpus[order(leaderspell, year)]

dt_colpus[, latent_personalism_lag := shift(latent_personalism, n = 1, type = "lag"), by = leaderspell]
dt_colpus[, rgdppc_lag := shift(rgdppc, n = 1, type = "lag"), by = leaderspell]
dt_colpus[, g_growth_rate_lag := shift(g_growth_rate, n = 1, type = "lag"), by = leaderspell]
dt_colpus[, oil_new_lag := shift(oil_new, n = 1, type = "lag"), by = leaderspell]
dt_colpus[, per_solider_lag := shift(per_solider, n = 1, type = "lag"), by = leaderspell]
dt_colpus[, mid_lag := shift(mid, n = 1, type = "lag"), by = leaderspell]
dt_colpus[, acd_intrastate_lag := shift(acd_intrastate, n = 1, type = "lag"), by = leaderspell]
dt_colpus[, exclpop_lag := shift(exclpop, n = 1, type = "lag"), by = leaderspell]
dt_colpus[, milethnic_hetero_lag := shift(milethnic_hetero, n = 1, type = "lag"), by = leaderspell]
dt_colpus[, defense_lag := shift(defense, n = 1, type = "lag"), by = leaderspell]
dt_colpus[, post_fail_lag := shift(post_fail_colpus, n = 1, type = "lag"), by = leaderspell]

min(dt_colpus$time_since_last_coup)
max(dt_colpus$time_since_last_coup)
mean(dt_colpus$time_since_last_coup)
sd(dt_colpus$time_since_last_coup)

################################ NAIVE OLS #####################################

naive_cb <- feols(counterbalancing ~ post_fail_colpus | leaderspell,
  cluster = "leaderspell", data = dt_colpus)
fixest::etable(naive_cb) #1%

naive_one <- feols(one_cw ~ post_fail_colpus | leaderspell,
  cluster = "leaderspell", data = dt_colpus)
fixest::etable(naive_one) #5%

naive_log <- feols(cb_log ~ post_fail_colpus | leaderspell,
  cluster = "leaderspell", data = dt_colpus)
fixest::etable(naive_log) #5%

################################## TWFE ########################################

#Control for Anticipation
twfe_cb <- feols(
  counterbalancing ~ post_fail_colpus + anticipation + t_leader +
    time_since_last_coup + latent_personalism_lag + rgdppc_lag + 
    g_growth_rate_lag + oil_new_lag + per_solider_lag + mid_lag + 
    acd_intrastate_lag + exclpop_lag + milethnic_hetero_lag + 
    defense_lag | leaderspell + year,
  cluster = "leaderspell", data = dt_colpus)
fixest::etable(twfe_cb) #5%

twfe_one <- feols(
  one_cw ~ post_fail_colpus + anticipation + t_leader +
    time_since_last_coup + latent_personalism_lag + rgdppc_lag + 
    g_growth_rate_lag + oil_new_lag + per_solider_lag + mid_lag + 
    acd_intrastate_lag + exclpop_lag + milethnic_hetero_lag + 
    defense_lag | leaderspell + year,
  cluster = "leaderspell", data = dt_colpus)
fixest::etable(twfe_one) #10%

twfe_log <- feols(
  cb_log ~ post_fail_colpus + anticipation + t_leader +
    time_since_last_coup + latent_personalism_lag + rgdppc_lag + 
    g_growth_rate_lag + oil_new_lag + per_solider_lag + mid_lag + 
    acd_intrastate_lag + exclpop_lag + milethnic_hetero_lag + 
    defense_lag | leaderspell + year,
  cluster = "leaderspell", data = dt_colpus)
fixest::etable(twfe_log) #10%

#Remove Anticipation (Test)
dt_colpus_ant <- dt_colpus %>%
  mutate(post_fail_colpus_ant = ifelse(duration == -1, 1, post_fail_colpus))

twfe_cb_ant <- feols(
  counterbalancing ~ post_fail_colpus + t_leader + time_since_last_coup +
    latent_personalism_lag + rgdppc_lag + g_growth_rate_lag + oil_new_lag +
    per_solider_lag + mid_lag + acd_intrastate_lag + exclpop_lag +
    milethnic_hetero_lag + defense_lag | leaderspell + year,
  cluster = "leaderspell", data = dt_colpus_ant)
fixest::etable(twfe_cb_ant) # 5%

twfe_one_ant <- feols(
  one_cw ~ post_fail_colpus + t_leader + time_since_last_coup +
    latent_personalism_lag + rgdppc_lag + g_growth_rate_lag + oil_new_lag +
    per_solider_lag + mid_lag + acd_intrastate_lag + exclpop_lag +
    milethnic_hetero_lag + defense_lag | leaderspell + year,
  cluster = "leaderspell", data = dt_colpus_ant)
fixest::etable(twfe_one_ant) # 5%

twfe_log_ant <- feols(
  cb_log ~ post_fail_colpus + t_leader + time_since_last_coup +
    latent_personalism_lag + rgdppc_lag + g_growth_rate_lag + oil_new_lag +
    per_solider_lag + mid_lag + acd_intrastate_lag + exclpop_lag +
    milethnic_hetero_lag + defense_lag | leaderspell + year,
  cluster = "leaderspell", data = dt_colpus_ant)
fixest::etable(twfe_log_ant) # 10%

################################### LDV TWFE ###################################

dt_colpus[, counterbalancing_lag := shift(counterbalancing, n = 1, type = "lag"), by = leaderspell]
dt_colpus[, one_cw_lag := shift(one_cw, n = 1, type = "lag"), by = leaderspell]
dt_colpus[, cb_log_lag := shift(cb_log, n = 1, type = "lag"), by = leaderspell]

twfe_cbl <- feols(
  counterbalancing ~ post_fail_colpus + t_leader + counterbalancing_lag + 
    time_since_last_coup + latent_personalism_lag + rgdppc_lag + 
    g_growth_rate_lag + oil_new_lag + per_solider_lag + mid_lag + 
    acd_intrastate_lag + exclpop_lag + milethnic_hetero_lag + 
    defense_lag | leaderspell + year,
  cluster = "leaderspell", data = dt_colpus)
fixest::etable(twfe_cbl) #5%

twfe_onel <- feols(
  one_cw ~ post_fail_colpus + t_leader + one_cw_lag + 
    time_since_last_coup + latent_personalism_lag + rgdppc_lag + 
    g_growth_rate_lag + oil_new_lag + per_solider_lag + mid_lag + 
    acd_intrastate_lag + exclpop_lag + milethnic_hetero_lag + 
    defense_lag | leaderspell + year,
  cluster = "leaderspell", data = dt_colpus)
fixest::etable(twfe_onel) #5%

twfe_ldv <- feols(
  cb_log ~ post_fail_colpus + t_leader + cb_log_lag + 
    time_since_last_coup + latent_personalism_lag + rgdppc_lag + 
    g_growth_rate_lag + oil_new_lag + per_solider_lag + mid_lag + 
    acd_intrastate_lag + exclpop_lag + milethnic_hetero_lag + 
    defense_lag | leaderspell + year,
  cluster = "leaderspell", data = dt_colpus)
fixest::etable(twfe_ldv) #10%

############################### Two-Stage DiD ##################################

did_cb <- did2s(
  data = dt_colpus,
  yname = "counterbalancing",
  treatment = "post_fail_colpus",
  first_stage = ~ t_leader + anticipation + time_since_last_coup + 
    latent_personalism_lag + rgdppc_lag + g_growth_rate_lag +
    oil_new_lag + per_solider_lag + mid_lag + acd_intrastate_lag +
    exclpop_lag + milethnic_hetero_lag + defense_lag | leaderspell + year,
  second_stage = ~ i(post_fail_colpus, ref = 0) +
    t_leader + anticipation + time_since_last_coup + latent_personalism_lag + 
    rgdppc_lag + g_growth_rate_lag + oil_new_lag + per_solider_lag + mid_lag + 
    acd_intrastate_lag + exclpop_lag + milethnic_hetero_lag + defense_lag,
  cluster_var = "leaderspell",
  bootstrap = FALSE,
  verbose = TRUE)
fixest::etable(did_cb) #5%

did_one <- did2s(
  data = dt_colpus,
  yname = "one_cw",
  treatment = "post_fail_colpus",
  first_stage = ~ t_leader + anticipation + time_since_last_coup + 
    latent_personalism_lag + rgdppc_lag + g_growth_rate_lag +
    oil_new_lag + per_solider_lag + mid_lag + acd_intrastate_lag +
    exclpop_lag + milethnic_hetero_lag + defense_lag | leaderspell + year,
  second_stage = ~ i(post_fail_colpus, ref = 0) +
    t_leader + anticipation + time_since_last_coup + latent_personalism_lag + 
    rgdppc_lag + g_growth_rate_lag + oil_new_lag + per_solider_lag + mid_lag + 
    acd_intrastate_lag + exclpop_lag + milethnic_hetero_lag + defense_lag,
  cluster_var = "leaderspell",
  bootstrap = FALSE,
  verbose = TRUE)
fixest::etable(did_one) #5%

did_log <- did2s(
  data = dt_colpus,
  yname = "cb_log",
  treatment = "post_fail_colpus",
  first_stage = ~ t_leader + anticipation + time_since_last_coup + 
    latent_personalism_lag + rgdppc_lag + g_growth_rate_lag +
    oil_new_lag + per_solider_lag + mid_lag + acd_intrastate_lag +
    exclpop_lag + milethnic_hetero_lag + defense_lag | leaderspell + year,
  second_stage = ~ i(post_fail_colpus, ref = 0) +
    t_leader + anticipation + time_since_last_coup + latent_personalism_lag + 
    rgdppc_lag + g_growth_rate_lag + oil_new_lag + per_solider_lag + mid_lag + 
    acd_intrastate_lag + exclpop_lag + milethnic_hetero_lag + defense_lag,
  cluster_var = "leaderspell",
  bootstrap = FALSE,
  verbose = TRUE)
fixest::etable(did_log) #5%
```

TWO-STAGE DID EVENT STUDY (FIGURE 5, FIGURE 6 AND FIGURE 7)

```{r}
############################### Counterbalancing ###############################

twfe_es <- feols(
  counterbalancing ~ i(duration, ref = c(-Inf, -1)) +
    t_leader + time_since_last_coup + latent_personalism_lag + rgdppc_lag + 
    g_growth_rate_lag + oil_new_lag + per_solider_lag + mid_lag +
    acd_intrastate_lag + exclpop_lag + milethnic_hetero_lag + defense_lag |
    leaderspell + year, data = dt_colpus)
fixest::etable(twfe_es)

#png("es_twfe_cb.png", width = 800, height = 500)
twfe_es |>
  fixest::iplot(
    main = "",
    xlab = "Time to Treatment",
    ylab = "Estimate and 95% CI",
    drop = c("[[:digit:]]{2}", "-Inf"),
    ref.line = -1,
    col = "#24496B")
axis(1, at = -9:9, labels = -9:9)
dev.off()

did2s_es <- did2s(
  data = dt_colpus,
  yname = "counterbalancing",
  treatment = "post_fail_colpus",
  first_stage = ~ t_leader + time_since_last_coup + latent_personalism_lag + 
    rgdppc_lag + g_growth_rate_lag + oil_new_lag + per_solider_lag + mid_lag + 
    acd_intrastate_lag + exclpop_lag + milethnic_hetero_lag + defense_lag | 
    leaderspell + year,
  second_stage = ~ i(duration, ref = c(-1, -Inf)) +
    t_leader + time_since_last_coup + latent_personalism_lag + rgdppc_lag +
    g_growth_rate_lag + oil_new_lag + per_solider_lag + mid_lag + 
    acd_intrastate_lag + exclpop_lag + milethnic_hetero_lag + defense_lag,
  cluster_var = "leaderspell")
fixest::etable(did2s_es)

#png("es_cb.png", width = 800, height = 500)
par(col.lab = "#24496B")
did2s_es |>
  fixest::iplot(
    main = "",
    xlab = "Time to Treatment",
    ylab = "Estimate and 95% CI",
    drop = c("[[:digit:]]{2}", "-Inf"),
    ref.line = -1,
    col = "#24496B",
    xlim = c(-9, 9))
axis(1, at = -9:9, labels = -9:9)
dev.off()

#Re-create in GGPLOT
iplot_cb <- fixest::iplot(did2s_es,
  plot = FALSE, drop = c("[[:digit:]]{2}", "-Inf"))$prms

ggplot(iplot_cb, aes(x = x, y = y)) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray") +
  geom_vline(xintercept = -1, lty = 2) + au_theme +
  geom_point(color = "#24496B", size = 1) +
  geom_errorbar(aes(ymin = ci_low, ymax = ci_high),
    width = 0.2, size = 0.4, color = "#24496B") +
  labs(x = "Time to Treatment",
       y = "Estimate and 95% CI") +
  scale_x_continuous(breaks = seq(-9, 9, 1)) +
  scale_y_continuous(breaks = seq(-0.5, 0.8, 0.1))
ggsave("did2sb_cb.png", width = 8, height = 4.5, dpi = 300)

############################## One Counterweight ###############################

twfe2_es <- feols(
  one_cw ~ i(duration, ref = c(-Inf, -1)) +
    t_leader + latent_personalism_lag + time_since_last_coup + rgdppc_lag + 
    g_growth_rate_lag + oil_new_lag + per_solider_lag + mid_lag +
    acd_intrastate_lag + exclpop_lag + milethnic_hetero_lag + defense_lag | 
    leaderspell + year, data = dt_colpus)
fixest::etable(twfe2_es)

#png("es_twfe_onecw.png", width = 800, height = 500)
par(col.lab = "#24496B")
twfe2_es |>
  fixest::iplot(
    main = "",
    xlab = "Time to Treatment",
    ylab = "Estimate and 95% CI",
    drop = c("[[:digit:]]{2}", "-Inf"),
    ref.line = -1,
    col = "#24496B")
axis(1, at = -9:9, labels = -9:9)
dev.off()

did2s2_es <- did2s(
  data = dt_colpus,
  yname = "one_cw",
  treatment = "post_fail_colpus",
  first_stage = ~ t_leader + time_since_last_coup + latent_personalism_lag + 
    rgdppc_lag + g_growth_rate_lag + oil_new_lag + per_solider_lag + mid_lag + 
    acd_intrastate_lag + exclpop_lag + milethnic_hetero_lag + defense_lag |
    leaderspell + year,
  second_stage = ~ i(duration, ref = c(-1, -Inf)) +
    t_leader + time_since_last_coup + latent_personalism_lag + rgdppc_lag + 
    g_growth_rate_lag + oil_new_lag + per_solider_lag + mid_lag + 
    acd_intrastate_lag + exclpop_lag + milethnic_hetero_lag + defense_lag,
  cluster_var = "leaderspell")
fixest::etable(did2s2_es)

png("es_onecw.png", width = 800, height = 500)
did2s2_es |>
  fixest::iplot(
    main = "",
    xlab = "Time to Treatment",
    ylab = "Estimate and 95% CI",
    drop = c("[[:digit:]]{2}", "-Inf"),
    ref.line = -1,
    col = "#24496B",
    ci_level = 0.95,
    xlim = c(-9, 9))
axis(1, at = -9:9, labels = -9:9)
dev.off()

# Re-create in GGPLOT
iplot_one <- fixest::iplot(did2s2_es,
  plot = FALSE, drop = c("[[:digit:]]{2}", "-Inf"))$prms

ggplot(iplot_one, aes(x = x, y = y)) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray") +
  geom_vline(xintercept = -1, lty = 2) + au_theme +
  geom_point(color = "#24496B", size = 1) +
  geom_errorbar(aes(ymin = ci_low, ymax = ci_high),
    width = 0.2, size = 0.4, color = "#24496B") +
  labs(x = "Time to Treatment",
       y = "Estimate and 95% CI") +
  scale_x_continuous(breaks = seq(-9, 9, 1)) +
  scale_y_continuous(breaks = seq(-0.5, 0.8, 0.1))
#ggsave("did2sb_one.png", width = 8, height = 4.5, dpi = 300)

############################## Log Counterweight ###############################

twfe3_es <- feols(
  cb_log ~ i(duration, ref = c(-Inf, -1)) +
    t_leader + time_since_last_coup + latent_personalism_lag + rgdppc_lag + 
    g_growth_rate_lag + oil_new_lag + per_solider_lag + mid_lag + 
    acd_intrastate_lag + exclpop_lag + milethnic_hetero_lag + defense_lag |
    leaderspell + year, data = dt_colpus)
fixest::etable(twfe3_es)

png("es_twfe_logcw.png", width = 800, height = 500)
par(col.lab = "#24496B")
twfe3_es |>
  fixest::iplot(
    main = "",
    xlab = "Time to Treatment",
    drop = c("[[:digit:]]{2}", "-Inf"),
    ref.line = -1,
    col = "#24496B")
axis(1, at = -9:9, labels = -9:9)
dev.off()

did2s3_es <- did2s(
  data = dt_colpus,
  yname = "cb_log",
  treatment = "post_fail_colpus",
  first_stage = ~ t_leader + time_since_last_coup + latent_personalism_lag + 
    rgdppc_lag + g_growth_rate_lag + oil_new_lag + per_solider_lag + mid_lag + 
    acd_intrastate_lag + exclpop_lag + milethnic_hetero_lag + defense_lag | 
    leaderspell + year,
  second_stage = ~ i(duration, ref = c(-1, -Inf)) +
    t_leader + time_since_last_coup + latent_personalism_lag + rgdppc_lag + 
    g_growth_rate_lag + oil_new_lag + per_solider_lag + mid_lag + 
    acd_intrastate_lag + exclpop_lag + milethnic_hetero_lag + defense_lag,
  cluster_var = "leaderspell")
fixest::etable(did2s3_es)

#png("es_logcw.png", width = 800, height = 500)
did2s3_es |>
  fixest::iplot(
    main = "",
    xlab = "Time to Treatment",
    ylab = "Estimate and 95% CI",
    drop = c("[[:digit:]]{2}", "-Inf"),
    ref.line = -1,
    col = "#24496B",
    ci_level = 0.95,
    xlim = c(-9, 9))
axis(1, at = -9:9, labels = -9:9)
dev.off()

# Re-create in GGPLOT
iplot_log <- fixest::iplot(did2s3_es,
  plot = FALSE,drop = c("[[:digit:]]{2}", "-Inf"))$prms

ggplot(iplot_log, aes(x = x, y = y)) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray") +
  geom_vline(xintercept = -1, lty = 2) + au_theme +
  geom_point(color = "#24496B", size = 1) +
  geom_errorbar(aes(ymin = ci_low, ymax = ci_high),
    width = 0.2, size = 0.4, color = "#24496B") +
  labs(x = "Time to Treatment",
       y = "Estimate and 95% CI") +
  scale_x_continuous(breaks = seq(-9, 9, 1)) +
  scale_y_continuous(breaks = seq(-0.5, 0.8, 0.1))
#ggsave("did2sb_log.png", width = 8, height = 4.5, dpi = 300)
```

APPENDIX: ALTERNATIVE DID (EVENT-STUDY) SPECIFICATIONS

```{r}
#Create Cohort Variable
did_colpus <- dt_colpus %>%
  mutate(leaderspell_id = as.numeric(as.factor(leaderspell)),
         year_id = as.numeric(as.character(year))) %>%
  arrange(leaderspell_id, year_id)

did_colpus <- did_colpus %>%
  group_by(leaderspell_id) %>%
  mutate(first_treat = if (any(post_fail_colpus == 1, na.rm = TRUE)) {
    min(year_id[post_fail_colpus == 1], na.rm = TRUE)
  } else {
    NA_integer_
  }) %>%
ungroup()

did_colpus$first_treat[is.na(did_colpus$first_treat)] <- 0

######################## STANDARD TWO-WAY FIXED EFFECTS ########################

#Log Counterweight
twfe_log2 <- event_study(
  data = did_colpus,
  yname = "cb_log",
  idname = "leaderspell_id",
  tname = "year_id",
  gname = "first_treat",
  xformla = ~ t_leader + time_since_last_coup + latent_personalism_lag +
    rgdppc_lag + g_growth_rate_lag + oil_new_lag + per_solider_lag + mid_lag +
    acd_intrastate_lag + exclpop_lag + milethnic_hetero_lag + defense_lag,
  estimator = "TWFE")

plot_event_study(twfe_log2, horizon = c(-9, 9))

twfe_log2 %>%
  filter(term %in% c(-9:9)) %>%
  ggplot(aes(x = term, y = estimate)) +
  geom_point(color = "#24496B", size = 1) +
  geom_errorbar(aes(ymin = estimate - 1.96 * std.error,
                    ymax = estimate + 1.96 * std.error),
    width = 0.2, size = 0.4, color = "#24496B") + au_theme +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray") +
  labs(x = "Time to Treatment",
       y = "Estimate and 95% CI") +
  scale_x_continuous(breaks = seq(-9, 9, 1)) +
  scale_y_continuous(breaks = seq(-0.5, 0.8, 0.1)) +
  geom_vline(xintercept = -1, lty = 2)
#ggsave("twfe2_log.png", width = 8, height = 4.5)

#One Counterweight
twfe_one2 <- event_study(
  data = did_colpus,
  yname = "one_cw",
  idname = "leaderspell_id",
  tname = "year_id",
  gname = "first_treat",
  xformla = ~ t_leader + time_since_last_coup + latent_personalism_lag +
    rgdppc_lag + g_growth_rate_lag + oil_new_lag + per_solider_lag + mid_lag +
    acd_intrastate_lag + exclpop_lag + milethnic_hetero_lag + defense_lag,
  estimator = "TWFE")

plot_event_study(twfe_one2, horizon = c(-9, 9))

twfe_one2 %>%
  filter(term %in% c(-9:9)) %>%
  ggplot(aes(x = term, y = estimate)) +
  geom_point(color = "#24496B", size = 1) +
  geom_errorbar(aes(ymin = estimate - 1.96 * std.error,
                    ymax = estimate + 1.96 * std.error),
    width = 0.2, size = 0.4, color = "#24496B") + au_theme +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray") +
  labs(x = "Time to Treatment",
       y = "Estimate and 95% CI") +
  scale_x_continuous(breaks = seq(-9, 9, 1)) +
  scale_y_continuous(breaks = seq(-0.5, 0.8, 0.1)) +
  geom_vline(xintercept = -1, lty = 2)
#ggsave("twfe2_one.png", width = 8, height = 4.5)

#Mean Counterbalancing
twfe_cb2 <- event_study(
  data = did_colpus,
  yname = "counterbalancing",
  idname = "leaderspell_id",
  tname = "year_id",
  gname = "first_treat",
  xformla = ~ t_leader + time_since_last_coup + latent_personalism_lag +
    rgdppc_lag + g_growth_rate_lag + oil_new_lag + per_solider_lag + mid_lag +
    acd_intrastate_lag + exclpop_lag + milethnic_hetero_lag + defense_lag,
  estimator = "TWFE")

plot_event_study(twfe_cb2, horizon = c(-9, 9))

twfe_cb2 %>%
  filter(term %in% c(-9:9)) %>%
  ggplot(aes(x = term, y = estimate)) +
  geom_point(color = "#24496B", size = 1) +
  geom_errorbar(aes(ymin = estimate - 1.96 * std.error,
                    ymax = estimate + 1.96 * std.error),
    width = 0.2, size = 0.4, color = "#24496B") + au_theme +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray") +
  labs(x = "Time to Treatment",
       y = "Estimate and 95% CI") +
  scale_x_continuous(breaks = seq(-9, 9, 1)) +
  scale_y_continuous(breaks = seq(-0.5, 0.8, 0.1)) +
  geom_vline(xintercept = -1, lty = 2)
#ggsave("twfe2_cb.png", width = 8, height = 4.5)

##################### BORUSYAK, JARAVEL AND SPIESS (2024) ######################

#Log Counterweight
bor_log <- event_study(
  data = did_colpus,
  yname = "cb_log",
  idname = "leaderspell_id",
  tname = "year_id",
  gname = "first_treat",
  xformla = ~ t_leader + time_since_last_coup + latent_personalism_lag +
    rgdppc_lag + g_growth_rate_lag + oil_new_lag + per_solider_lag + mid_lag +
    acd_intrastate_lag + exclpop_lag + milethnic_hetero_lag + defense_lag,
  estimator = "impute")

plot_event_study(bor_log, horizon = c(-9, 9))

bor_log %>%
  filter(term %in% c(-9:9)) %>%
  ggplot(aes(x = term, y = estimate)) +
  geom_point(color = "#24496B", size = 1) +
  geom_errorbar(aes(ymin = estimate - 1.96 * std.error,
                    ymax = estimate + 1.96 * std.error),
    width = 0.2, size = 0.4, color = "#24496B") + au_theme +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray") +
  labs(x = "Time to Treatment", y = "Estimate and 95% CI") +
  scale_x_continuous(breaks = seq(-9, 9, 1)) +
  scale_y_continuous(breaks = seq(-0.5, 0.8, 0.1)) +
  geom_vline(xintercept = -1, lty = 2)
#ggsave("bor2_log.png", width = 8, height = 4.5)

#One Counterweight
bor_one <- event_study(
  data = did_colpus,
  yname = "one_cw",
  idname = "leaderspell_id",
  tname = "year_id",
  gname = "first_treat",
  xformla = ~ t_leader + time_since_last_coup + latent_personalism_lag +
    rgdppc_lag + g_growth_rate_lag + oil_new_lag + per_solider_lag + mid_lag +
    acd_intrastate_lag + exclpop_lag + milethnic_hetero_lag + defense_lag,
  estimator = "impute")

plot_event_study(bor_one, horizon = c(-9, 9))

bor_one %>%
  filter(term %in% c(-9:9)) %>%
  ggplot(aes(x = term, y = estimate)) +
  geom_point(color = "#24496B", size = 1) +
  geom_errorbar(aes(ymin = estimate - 1.96 * std.error,
                    ymax = estimate + 1.96 * std.error),
    width = 0.2, size = 0.4, color = "#24496B") + au_theme +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray") +
  labs( x = "Time to Treatment",
        y = "Estimate and 95% CI") +
  scale_x_continuous(breaks = seq(-9, 9, 1)) +
  scale_y_continuous(breaks = seq(-0.5, 0.8, 0.1)) +
  geom_vline(xintercept = -1, lty = 2)
#ggsave("bor2_one.png", width = 8, height = 4.5)

#Mean Counterbalancing
bor_cb <- event_study(
  data = did_colpus,
  yname = "counterbalancing",
  idname = "leaderspell_id",
  tname = "year_id",
  gname = "first_treat",
  xformla = ~ t_leader + time_since_last_coup + latent_personalism_lag +
    rgdppc_lag + g_growth_rate_lag + oil_new_lag + per_solider_lag + mid_lag +
    acd_intrastate_lag + exclpop_lag + milethnic_hetero_lag + defense_lag,
  estimator = "impute")

plot_event_study(bor_cb, horizon = c(-9, 9))

bor_cb %>%
  filter(term %in% c(-9:9)) %>%
  ggplot(aes(x = term, y = estimate)) +
  geom_point(color = "#24496B", size = 1) +
  geom_errorbar(aes(ymin = estimate - 1.96 * std.error,
                    ymax = estimate + 1.96 * std.error),
    width = 0.2, size = 0.4, color = "#24496B") + au_theme +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray") +
  labs(x = "Time to Treatment",
       y = "Estimate and 95% CI") +
  scale_x_continuous(breaks = seq(-9, 9, 1)) +
  scale_y_continuous(breaks = seq(-0.5, 0.8, 0.1)) +
  geom_vline(xintercept = -1, lty = 2)
#ggsave("bor2_cb.png", width = 8, height = 4.5)

######################## GARDNER (2021) FIRST CONTROL ##########################

#Log Counterweight
did2s2_log <- event_study(
  data = did_colpus,
  yname = "cb_log",
  idname = "leaderspell_id",
  tname = "year_id",
  gname = "first_treat",
  xformla = ~ t_leader + time_since_last_coup + latent_personalism_lag +
    rgdppc_lag + g_growth_rate_lag + oil_new_lag + per_solider_lag + mid_lag +
    acd_intrastate_lag + exclpop_lag + milethnic_hetero_lag + defense_lag,
  estimator = "did2s")

plot_event_study(did2s2_log, horizon = c(-9, 9))

did2s2_log %>%
  filter(term %in% c(-9:9) & term != -1) %>%
  ggplot(aes(x = term, y = estimate)) +
  geom_point(color = "#24496B", size = 1) +
  geom_errorbar(aes(ymin = estimate - 1.96 * std.error,
                    ymax = estimate + 1.96 * std.error),
    width = 0.2, size = 0.4, color = "#24496B") + au_theme +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray") +
  labs(x = "Time to Treatment",
       y = "Estimate and 95% CI") +
  scale_x_continuous(breaks = seq(-9, 9, 1)) +
  scale_y_continuous(breaks = seq(-0.5, 0.8, 0.1)) +
  geom_vline(xintercept = -1, lty = 2)
#ggsave("did2s2_log.png", width = 8, height = 4.5)

#One Counterweight
did2s2_one <- event_study(
  data = did_colpus,
  yname = "one_cw",
  idname = "leaderspell_id",
  tname = "year_id",
  gname = "first_treat",
  xformla = ~ t_leader + time_since_last_coup + latent_personalism_lag +
    rgdppc_lag + g_growth_rate_lag + oil_new_lag + per_solider_lag + mid_lag +
    acd_intrastate_lag + exclpop_lag + milethnic_hetero_lag + defense_lag,
  estimator = "did2s")

plot_event_study(did2s2_one, horizon = c(-9, 9))

did2s2_one %>%
  filter(term %in% c(-9:9) & term != -1) %>%
  ggplot(aes(x = term, y = estimate)) +
  geom_point(color = "#24496B", size = 1) +
  geom_errorbar(aes(ymin = estimate - 1.96 * std.error,
                    ymax = estimate + 1.96 * std.error),
    width = 0.2, size = 0.4, color = "#24496B") + au_theme +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray") +
  labs(x = "Time to Treatment",
       y = "Estimate and 95% CI") +
  scale_x_continuous(breaks = seq(-9, 9, 1)) +
  scale_y_continuous(breaks = seq(-0.5, 0.8, 0.1)) +
  geom_vline(xintercept = -1, lty = 2)
#ggsave("did2s2_one.png", width = 8, height = 4.5)

# Mean Counterbalancing
did2s2_cb <- event_study(
  data = did_colpus,
  yname = "counterbalancing",
  idname = "leaderspell_id",
  tname = "year_id",
  gname = "first_treat",
  xformla = ~ t_leader + time_since_last_coup + latent_personalism_lag +
    rgdppc_lag + g_growth_rate_lag + oil_new_lag + per_solider_lag + mid_lag +
    acd_intrastate_lag + exclpop_lag + milethnic_hetero_lag + defense_lag,
  estimator = "did2s")

plot_event_study(did2s2_cb, horizon = c(-9, 9))

did2s2_cb %>%
  filter(term %in% c(-9:9) & term != -1) %>%
  ggplot(aes(x = term, y = estimate)) +
  geom_point(color = "#24496B", size = 1) +
  geom_errorbar(aes(ymin = estimate - 1.96 * std.error,
                    ymax = estimate + 1.96 * std.error),
    width = 0.2, size = 0.4, color = "#24496B") + au_theme +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray") +
  labs(x = "Time to Treatment",
       y = "Estimate and 95% CI") +
  scale_x_continuous(breaks = seq(-9, 9, 1)) +
  scale_y_continuous(breaks = seq(-0.5, 0.8, 0.1)) +
  geom_vline(xintercept = -1, lty = 2)
#ggsave("did2s2_cb.png", width = 8, height = 4.5)

######################## Callaway and Sant'Anna (2020) #########################

#Log Counterweights
att_gt_log <- att_gt(
  yname = "cb_log",
  tname = "year_id",
  idname = "leaderspell_id",
  gname = "first_treat",
  data = did_colpus,
  est_method = "reg",
  panel = FALSE,
  control_group = "notyettreated",
  clustervars = "leaderspell_id")

es_att_log <- aggte(att_gt_log, type = "dynamic", na.rm = T)

df_att_log <- data.frame(
  event_time = es_att_log$egt,
  att = es_att_log$att.egt,
  ci_lower = es_att_log$att.egt - 1.96 * es_att_log$se.egt,
  ci_upper = es_att_log$att.egt + 1.96 * es_att_log$se.egt)

df_att_log %>%
  filter(event_time %in% c(-9:9) & event_time != -1) %>%
  ggplot(aes(x = event_time, y = att)) +
  geom_point(color = "#24496B", size = 1) +
  geom_errorbar(aes(ymin = ci_lower, ymax = ci_upper),
    width = 0.2, size = 0.4, color = "#24496B") + au_theme +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray") +
  labs(x = "Time to Treatment",
       y = "Estimate and 95% CI") +
  scale_x_continuous(breaks = seq(-9, 9, 1)) +
  scale_y_continuous(breaks = seq(-0.5, 0.8, 0.1)) +
  geom_vline(xintercept = -1, lty = 2)
#ggsave("gt_log2.png", width = 8, height = 4.5)

#One Counterweight
att_gt_one <- att_gt(
  yname = "one_cw",
  tname = "year_id",
  idname = "leaderspell_id",
  gname = "first_treat",
  data = did_colpus,
  est_method = "reg",
  panel = FALSE,
  control_group = "notyettreated",
  clustervars = "leaderspell_id")

es_att_one <- aggte(att_gt_one, type = "dynamic", na.rm = T)

df_att_one <- data.frame(
  event_time = es_att_one$egt,
  att = es_att_one$att.egt,
  ci_lower = es_att_one$att.egt - 1.96 * es_att_one$se.egt,
  ci_upper = es_att_one$att.egt + 1.96 * es_att_one$se.egt)

df_att_one %>%
  filter(event_time %in% c(-9:9) & event_time != -1) %>%
  ggplot(aes(x = event_time, y = att)) +
  geom_point(color = "#24496B", size = 1) +
  geom_errorbar(aes(ymin = ci_lower,ymax = ci_upper),
    width = 0.2, size = 0.4, color = "#24496B") + au_theme +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray") +
  labs(x = "Time to Treatment",
       y = "Estimate and 95% CI") +
  scale_x_continuous(breaks = seq(-9, 9, 1)) +
  scale_y_continuous(breaks = seq(-0.5, 0.8, 0.1)) +
  geom_vline(xintercept = -1, lty = 2)
#ggsave("gt_one2.png", width = 8, height = 4.5)

# Mean Counterbalancing
att_gt_cb <- att_gt(
  yname = "counterbalancing",
  tname = "year_id",
  idname = "leaderspell_id",
  gname = "first_treat",
  data = did_colpus,
  est_method = "reg",
  panel = FALSE,
  control_group = "notyettreated",
  clustervars = "leaderspell_id")

es_att_cb <- aggte(att_gt_cb, type = "dynamic", na.rm = T)

df_att_cb <- data.frame(
  event_time = es_att_cb$egt,
  att = es_att_cb$att.egt,
  ci_lower = es_att_cb$att.egt - 1.96 * es_att_cb$se.egt,
  ci_upper = es_att_cb$att.egt + 1.96 * es_att_cb$se.egt)

df_att_cb %>%
  filter(event_time %in% c(-9:9) & event_time != -1) %>%
  ggplot(aes(x = event_time, y = att)) +
  geom_point(color = "#24496B", size = 1) +
  geom_errorbar(aes(ymin = ci_lower,ymax = ci_upper),
    width = 0.2, size = 0.4, color = "#24496B") + au_theme +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray") +
  labs(x = "Time to Treatment",
       y = "Estimate and 95% CI") +
  scale_x_continuous(breaks = seq(-9, 9, 1)) +
  scale_y_continuous(breaks = seq(-0.5, 0.8, 0.1)) +
  geom_vline(xintercept = -1, lty = 2)
#ggsave("gt_cb2.png", width = 8, height = 4.5)
```

APPENDIX: COEFFICIENT PLOT COUNTERBALANCING (TABLE A4.1)

```{r}
#Naive (95%)
coef_naive <- naive_cb$coefficients[1]
se_naive <- naive_cb$se[1]
ci_naive_low <- coef_naive - (1.96 * se_naive)
ci_naive_high <- coef_naive + (1.96 * se_naive)

#TWFE (95%)
coef_twfe <- twfe_cb$coefficients[1]
se_twfe <- twfe_cb$se[1]
ci_twfe_low <- coef_twfe - (1.96 * se_twfe)
ci_twfe_high <- coef_twfe + (1.96 * se_twfe)

#Two-Stage Did (95%)
coef_ts <- did_cb$coefficients[1]
se_ts <- did_cb$se[1]
ci_ts_low <- coef_ts - (1.96 * se_ts)
ci_ts_high <- coef_ts + (1.96 * se_ts)

#Naive (90%)
ci_naive_low90 <- coef_naive - (1.645 * se_naive)
ci_naive_high90 <- coef_naive + (1.645 * se_naive)

#TWFE (90%)
ci_twfe_low90 <- coef_twfe - (1.645 * se_twfe)
ci_twfe_high90 <- coef_twfe + (1.645 * se_twfe)

#Two-Stage Did (90%)
ci_ts_low90 <- coef_ts - (1.645 * se_ts)
ci_ts_high90 <- coef_ts + (1.645 * se_ts)

#Merge
row_labs <- c("mean", "low", "high", "low90", "high90")
naive_merge <- as.data.frame(cbind(
  coef_naive, ci_naive_low, ci_naive_high,
  ci_naive_low90, ci_naive_high90))
colnames(naive_merge) <- row_labs

twfe_merge <- as.data.frame(cbind(
  coef_twfe, ci_twfe_low, ci_twfe_high,
  ci_twfe_low90, ci_twfe_high90))
colnames(twfe_merge) <- row_labs

ts_merge <- as.data.frame(cbind(
  coef_ts, ci_ts_low, ci_ts_high,
  ci_ts_low90, ci_ts_high90))
colnames(ts_merge) <- row_labs

merge95 <- as.data.frame(rbind(naive_merge, twfe_merge, ts_merge))
merge95$type <- c("1Naïve", "2TWFE", "3Two-Stage DiD")
merge95$order <- 3:1

merge95$mean <- merge95$mean * 100
merge95$low <- merge95$low * 100
merge95$high <- merge95$high * 100
merge95$low90 <- merge95$low90 * 100
merge95$high90 <- merge95$high90 * 100
merge951 <- merge95

coef_labs <- c("Naïve", "TWFE", "2SDID")

ggplot(merge951, aes(y = order, shape = type)) +
  geom_pointrange(aes(x = mean, xmin = low, xmax = high), alpha = 0.25) +
  geom_pointrange(aes(x = mean, xmin = low90, xmax = high90)) +
  au_theme + geom_vline(xintercept = 0, lty = 2, color = "#838383") +
  # scale_x_continuous(breaks = seq(0, 40, 5)) +
  scale_shape_manual(values = merge951$order, labels = coef_labs) +
  labs(x = "Average Treatment Effect",
  caption = "Note: Gray Line Depicts 95% Confidence Intervals;
       Black Line Depicts 90% Confidence Intervals") +
  theme(axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        legend.title = element_blank(),
        plot.caption = element_text(size = 30, hjust = 0))
#ggsave("coefpolot_cb3.png", width = 8, height = 4.5)
```

APPENDIX: COEFFICIENT PLOT ONE COUNTERWEIGHT (TABLE A4.1)

```{r}
# Naive (95%)
coef_naive <- naive_one$coefficients[1]
se_naive <- naive_one$se[1]
ci_naive_low <- coef_naive - (1.96 * se_naive)
ci_naive_high <- coef_naive + (1.96 * se_naive)

# TWFE (95%)
coef_twfe <- twfe_one$coefficients[1]
se_twfe <- twfe_one$se[1]
ci_twfe_low <- coef_twfe - (1.96 * se_twfe)
ci_twfe_high <- coef_twfe + (1.96 * se_twfe)

# Two-Stage Did (95%)
coef_ts <- did_one$coefficients[1]
se_ts <- did_one$se[1]
ci_ts_low <- coef_ts - (1.96 * se_ts)
ci_ts_high <- coef_ts + (1.96 * se_ts)

# Naive (90%)
ci_naive_low90 <- coef_naive - (1.645 * se_naive)
ci_naive_high90 <- coef_naive + (1.645 * se_naive)

# TWFE (90%)
ci_twfe_low90 <- coef_twfe - (1.645 * se_twfe)
ci_twfe_high90 <- coef_twfe + (1.645 * se_twfe)

# Two-Stage Did (90%)
ci_ts_low90 <- coef_ts - (1.645 * se_ts)
ci_ts_high90 <- coef_ts + (1.645 * se_ts)

# Merge
row_labs <- c("mean", "low", "high", "low90", "high90")
naive_merge <- as.data.frame(cbind(
  coef_naive, ci_naive_low, ci_naive_high,
  ci_naive_low90, ci_naive_high90))
colnames(naive_merge) <- row_labs

twfe_merge <- as.data.frame(cbind(
  coef_twfe, ci_twfe_low, ci_twfe_high,
  ci_twfe_low90, ci_twfe_high90))
colnames(twfe_merge) <- row_labs

ts_merge <- as.data.frame(cbind(
  coef_ts, ci_ts_low, ci_ts_high,
  ci_ts_low90, ci_ts_high90))
colnames(ts_merge) <- row_labs

merge95 <- as.data.frame(rbind(naive_merge, twfe_merge, ts_merge))
merge95$type <- c("1Naïve", "2TWFE", "3Two-Stage DiD")
merge95$order <- 3:1

merge95$mean <- merge95$mean * 100
merge95$low <- merge95$low * 100
merge95$high <- merge95$high * 100
merge95$low90 <- merge95$low90 * 100
merge95$high90 <- merge95$high90 * 100
merge952 <- merge95

ggplot(merge952, aes(y = order, shape = type)) +
  geom_pointrange(aes(x = mean, xmin = low, xmax = high), alpha = 0.25) +
  geom_pointrange(aes(x = mean, xmin = low90, xmax = high90)) +
  au_theme + geom_vline(xintercept = 0, lty = 2, color = "#838383") +
  scale_x_continuous(breaks = seq(0, 40, 5)) +
  scale_shape_manual(values = merge952$order, labels = coef_labs) +
  labs(x = "Average Treatment Effect",
  caption = "Note: Gray Line Depicts 95% Confidence Intervals;
       Black Line Depicts 90% Confidence Intervals") +
  theme(axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        legend.title = element_blank(),
        plot.caption = element_text(size = 30, hjust = 0))
#ggsave("coefpolot_one3.png", width = 8, height = 4.5)
```

APPENDIX: COEFFICIENT PLOT LOG COUNTERWEIGHTS (TABLE A4.1) - EXPONENTIATED

```{r}
# Naive (95%)
coef_naive <- (exp(naive_log$coefficients[1]) - 1) # 0.1434
se_naive <- naive_log$se[1]
ci_naive_low <- coef_naive - (1.96 * se_naive)
ci_naive_high <- coef_naive + (1.96 * se_naive)

# TWFE (95%)
coef_twfe <- (exp(twfe_log$coefficients[1]) - 1) # 0.1453
se_twfe <- twfe_log$se[1]
ci_twfe_low <- coef_twfe - (1.96 * se_twfe)
ci_twfe_high <- coef_twfe + (1.96 * se_twfe)

# Two-Stage Did (95%)
coef_ts <- (exp(did_log$coefficients[1]) - 1) # 0.2577
se_ts <- did_log$se[1]
ci_ts_low <- coef_ts - (1.96 * se_ts)
ci_ts_high <- coef_ts + (1.96 * se_ts)

# Naive (90%)
ci_naive_low90 <- coef_naive - (1.645 * se_naive)
ci_naive_high90 <- coef_naive + (1.645 * se_naive)

# TWFE (90%)
ci_twfe_low90 <- coef_twfe - (1.645 * se_twfe)
ci_twfe_high90 <- coef_twfe + (1.645 * se_twfe)

# Two-Stage Did (90%)
ci_ts_low90 <- coef_ts - (1.645 * se_ts)
ci_ts_high90 <- coef_ts + (1.645 * se_ts)

# Merge
row_labs <- c("mean", "low", "high", "low90", "high90")
naive_merge <- as.data.frame(cbind(
  coef_naive, ci_naive_low, ci_naive_high,
  ci_naive_low90, ci_naive_high90))
colnames(naive_merge) <- row_labs

twfe_merge <- as.data.frame(cbind(
  coef_twfe, ci_twfe_low, ci_twfe_high,
  ci_twfe_low90, ci_twfe_high90))
colnames(twfe_merge) <- row_labs

ts_merge <- as.data.frame(cbind(
  coef_ts, ci_ts_low, ci_ts_high,
  ci_ts_low90, ci_ts_high90))
colnames(ts_merge) <- row_labs

merge95 <- as.data.frame(rbind(naive_merge, twfe_merge, ts_merge))
merge95$type <- c("1Naïve", "2TWFE", "3Two-Stage DiD")
merge95$order <- 3:1

merge95$mean <- merge95$mean * 100
merge95$low <- merge95$low * 100
merge95$high <- merge95$high * 100
merge95$low90 <- merge95$low90 * 100
merge95$high90 <- merge95$high90 * 100
merge953 <- merge95

coef_labs_new <- c("Naïve", "TWFE", "Two-Stage DiD")

ggplot(merge953, aes(y = order, shape = type)) +
  geom_pointrange(aes(x = mean, xmin = low, xmax = high), alpha = 0.25) +
  geom_pointrange(aes(x = mean, xmin = low90, xmax = high90)) +
  au_theme + geom_vline(xintercept = 0, lty = 2, color = "#838383") +
  # scale_x_continuous(breaks = seq(0, 40, 5)) +
  scale_shape_manual(values = merge953$order, labels = coef_labs_new) +
  labs(x = "Average Marginal Effect (%)",
  caption = "Note: Gray Line Depicts 95% Confidence Intervals;
  Black Line Depicts 90% Confidence Intervals") +
  theme(axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        legend.title = element_blank(),
        plot.caption = element_text(size = 30, hjust = 0))
#ggsave("coefpolot_log3.png", width = 8, height = 4.5)
```

APPENDIX: JOINT COEFFICIENT PLOT (FIGURE A4.1)

```{r}
# Facet Wrap Plot
merge951$outcome <- "3Counterbalancing"
merge952$outcome <- "2One Counterweight"
merge953$outcome <- "1Log Counterweights"

merge_joint <- as.data.frame(rbind(merge951, merge952, merge953))

ggplot(merge_joint, aes(y = order, shape = type)) +
  geom_pointrange(aes(x = mean, xmin = low, xmax = high), alpha = 0.25) +
  geom_pointrange(aes(x = mean, xmin = low90, xmax = high90)) +
  facet_wrap(~outcome, labeller = as_labeller(c(
    "3Counterbalancing" = "Mean Counterbalancing",
    "2One Counterweight" = "One Counterweight",
    "1Log Counterweights" = "Log Counterweights"))) +
  au_theme +  geom_vline(xintercept = 0, lty = 2, color = "#838383") +
  scale_shape_manual(values = merge_joint$order, labels = coef_labs) +
  # scale_x_continuous(breaks = seq(0, 40, 5)) +
  labs(x = "Estimate (ATET)",
    caption = "Note: Gray Line Depicts 95% Confidence Intervals; 
    Black Line Depicts 90% Confidence Intervals; 
    Exponentiated Coefficients Shown") +
  theme(axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        legend.title = element_blank(),
        legend.spacing = unit(1, "cm"),
        plot.caption = element_text(size = 28, hjust = 0),
        strip.background = element_rect(color = "white", fill = "white"),
        strip.text = element_text(size = 6, color = "#24496B"))
#ggsave("coefpolot_joint3.png", width = 8, height = 4.5)
```

APPENDIX: ROBUSTNESS CHECKS (ALTERNATIVE SPECIFICATIONS)

```{r}
############################## Logistic Regression #############################

# Naive (Counterbalancing)
dt_colpus_r <- dt_colpus %>% rename(cb = counterbalancing)

naive_new_glm1 <- feglm(cb ~ post_fail_colpus | leaderspell,
  cluster = "leaderspell", family = "logit", data = dt_colpus_r)
fixest::etable(naive_new_glm1)

# TWFE (Counterbalancing)
twfe_logit1 <- feglm(
  cb ~ post_fail_colpus + # anticipation +
    t_leader + time_since_last_coup + latent_personalism_lag +
    rgdppc_lag + g_growth_rate_lag + oil_new_lag + per_solider_lag +
    mid_lag + acd_intrastate_lag + exclpop_lag + milethnic_hetero_lag +
    defense_lag | leaderspell + year,
  cluster = "leaderspell",
  family = "logit", data = dt_colpus_r)
fixest::etable(twfe_logit1)

# Naive (One Counterweight)
naive_new_glm2 <- feglm(one_cw ~ post_fail_colpus | leaderspell,
  cluster = "leaderspell", family = "logit", data = dt_colpus_r)
fixest::etable(naive_new_glm2)

# TWFE (One Counterweight)
twfe_logit2 <- feglm(
  one_cw ~ post_fail_colpus + # anticipation +
    t_leader + time_since_last_coup + latent_personalism_lag +
    rgdppc_lag + g_growth_rate_lag + oil_new_lag + per_solider_lag +
    mid_lag + acd_intrastate_lag + exclpop_lag + milethnic_hetero_lag +
    defense_lag | leaderspell + year,
  cluster = "leaderspell",
  family = "logit", data = dt_colpus_r)
fixest::etable(twfe_logit2)

etable_output <- fixest::etable(
  list(naive_new_glm1, twfe_logit1, naive_new_glm2, twfe_logit2),
  headers = c("Naïve", "TWFE", "Naïve", "TWFE"))

custom_css <- "<style>
table.classic-table {
  border-collapse: collapse;
  margin: auto;
  font-size: 10px;
  font-family: 'Arial', sans-serif;
  width: auto;
}

table.classic-table th,
table.classic-table td {
  padding: 4px 6px;
  text-align: center;
  border-bottom: 1px solid #ccc;
}

table.classic-table th {
  border-bottom: 2px solid #444;
  font-weight: bold;
  background-color: #f9f9f9;
}

table.classic-table caption {
  caption-side: top;
  font-weight: bold;
  font-size: 12px;
  margin-bottom: 6px;
  color: #222;
}
</style>"

kable_log <- kable(
  etable_output,
  format = "html",
  caption = "",
  escape = FALSE
) %>%
  kable_styling(
    full_width = FALSE,
    position = "center",
    font_size = 14,
    htmltable_class = "classic-table"
  ) %>%
  footnote(
    general = "Significance levels: *** p < 0.001, ** p < 0.01, * p < 0.05, . p < 0.1.",
    general_title = "Note:")

full_html <- paste0(custom_css, as.character(kable_log))
#writeLines(full_html, "kable_log.html")

############################## GWF Regime Control ##############################

# GWF Regime as Control
twfe_reg <- (feols(
  cb ~ post_fail_colpus + anticipation + t_leader + time_since_last_coup +
    gwf_personal + rgdppc_lag + g_growth_rate_lag + oil_new_lag +
    per_solider_lag + mid_lag + acd_intrastate_lag + exclpop_lag +
    milethnic_hetero_lag + defense_lag | leaderspell + year,
  cluster = "leaderspell", data = dt_colpus_r))
fixest::etable(twfe_reg)

twfe_reg2 <- (feols(
  one_cw ~ post_fail_colpus + anticipation + t_leader + time_since_last_coup +
    gwf_military + gwf_party + rgdppc_lag + g_growth_rate_lag + oil_new_lag +
    per_solider_lag + mid_lag + acd_intrastate_lag + exclpop_lag +
    milethnic_hetero_lag + defense_lag | leaderspell + year,
  cluster = "leaderspell", data = dt_colpus_r))
fixest::etable(twfe_reg2)

twfe_reg3 <- (feols(
  cb_log ~ post_fail_colpus + anticipation + t_leader + time_since_last_coup +
    gwf_military + gwf_party + rgdppc_lag + g_growth_rate_lag + oil_new_lag +
    per_solider_lag + mid_lag + acd_intrastate_lag + exclpop_lag +
    milethnic_hetero_lag + defense_lag | leaderspell + year,
  cluster = "leaderspell", data = dt_colpus_r))
fixest::etable(twfe_reg3)

etable_output <- fixest::etable(
  list(twfe_reg, twfe_reg2, twfe_reg3),
  headers = c("TWFE", "TWFE", "TWFE"))

kable_g <- kable(
  etable_output,
  format = "html",
  caption = "",
  escape = FALSE
) %>%
  kable_styling(
    full_width = FALSE,
    position = "center",
    font_size = 14,
    htmltable_class = "classic-table"
  ) %>%
  footnote(
    general = "Significance levels: *** p < 0.001, ** p < 0.01, * p < 0.05, . p < 0.1.",
    general_title = "Note:")

full_html <- paste0(custom_css, as.character(kable_g))
#writeLines(full_html, "kable_g.html")
#Redundant if I control for leader fixed-effects

############################ Cluster Standard Errors ###########################

#Regime Clustered Standard Error
naive_cb_r <- feols(cb ~ post_fail_colpus | leaderspell,
  cluster = "gwf_casename", data = dt_colpus_r) # 5%
fixest::etable(naive_cb_r)

naive_one_r <- feols(one_cw ~ post_fail_colpus | leaderspell,
  cluster = "gwf_casename", data = dt_colpus_r) # 5%
fixest::etable(naive_one_r)

naive_log_r <- feols(cb_log ~ post_fail_colpus | leaderspell,
  cluster = "gwf_casename", data = dt_colpus_r) # 5%
fixest::etable(naive_log_r)

twfe_cb_r <- feols(
  cb ~ post_fail_colpus + anticipation +
    t_leader + time_since_last_coup + latent_personalism_lag +
    rgdppc_lag + g_growth_rate_lag + oil_new_lag + per_solider_lag +
    mid_lag + acd_intrastate_lag + exclpop_lag + milethnic_hetero_lag +
    defense_lag | leaderspell + year,
  cluster = "gwf_casename",
  data = dt_colpus_r) # 5%
fixest::etable(twfe_cb_r)

twfe_one_r <- feols(
  one_cw ~ post_fail_colpus + anticipation +
    t_leader + time_since_last_coup + latent_personalism_lag +
    rgdppc_lag + g_growth_rate_lag + oil_new_lag + per_solider_lag +
    mid_lag + acd_intrastate_lag + exclpop_lag + milethnic_hetero_lag +
    defense_lag | leaderspell + year,
  cluster = "gwf_casename",
  data = dt_colpus_r) # 5%
fixest::etable(twfe_one_r)

twfe_log_r <- feols(
  cb_log ~ post_fail_colpus + anticipation +
    t_leader + time_since_last_coup + latent_personalism_lag +
    rgdppc_lag + g_growth_rate_lag + oil_new_lag + per_solider_lag +
    mid_lag + acd_intrastate_lag + exclpop_lag + milethnic_hetero_lag +
    defense_lag | leaderspell + year,
  cluster = "gwf_casename",
  data = dt_colpus_r) # 5%
fixest::etable(twfe_log_r)

did_cb_r <- did2s(
  data = dt_colpus_r,
  yname = "cb",
  treatment = "post_fail_colpus",
  first_stage = ~ t_leader + anticipation + latent_personalism_lag + 
    rgdppc_lag + g_growth_rate_lag + oil_new_lag + per_solider_lag + 
    mid_lag + acd_intrastate_lag + exclpop_lag + milethnic_hetero_lag + 
    defense_lag | leaderspell + year,
  second_stage = ~ i(post_fail_colpus, ref = 0) +
    t_leader + anticipation + latent_personalism_lag + rgdppc_lag + 
    g_growth_rate_lag + oil_new_lag + per_solider_lag + mid_lag + 
    acd_intrastate_lag + exclpop_lag + milethnic_hetero_lag + defense_lag,
  cluster_var = "gwf_casename",
  bootstrap = FALSE,
  verbose = TRUE) # 5%
fixest::etable(did_cb_r)

did_one_r <- did2s(
  data = dt_colpus_r,
  yname = "one_cw",
  treatment = "post_fail_colpus",
  first_stage = ~ t_leader + anticipation + latent_personalism_lag + 
    rgdppc_lag + g_growth_rate_lag + oil_new_lag + per_solider_lag + 
    mid_lag + acd_intrastate_lag + exclpop_lag + milethnic_hetero_lag + 
    defense_lag | leaderspell + year,
  second_stage = ~ i(post_fail_colpus, ref = 0) +
    t_leader + anticipation + latent_personalism_lag + rgdppc_lag + 
    g_growth_rate_lag + oil_new_lag + per_solider_lag + mid_lag + 
    acd_intrastate_lag + exclpop_lag + milethnic_hetero_lag + defense_lag,
  cluster_var = "gwf_casename",
  bootstrap = FALSE,
  verbose = TRUE) # 5%
fixest::etable(did_one_r)

did_log_r <- did2s(
  data = dt_colpus_r,
  yname = "cb_log",
  treatment = "post_fail_colpus",
  first_stage = ~ t_leader + anticipation + latent_personalism_lag + 
    rgdppc_lag + g_growth_rate_lag + oil_new_lag + per_solider_lag +
    mid_lag + acd_intrastate_lag + exclpop_lag + milethnic_hetero_lag +
    defense_lag | leaderspell + year,
  second_stage = ~ i(post_fail_colpus, ref = 0) +
    t_leader + anticipation + latent_personalism_lag + rgdppc_lag +
    g_growth_rate_lag + oil_new_lag + per_solider_lag + mid_lag +
    acd_intrastate_lag + exclpop_lag + milethnic_hetero_lag + defense_lag,
  cluster_var = "gwf_casename",
  bootstrap = FALSE,
  verbose = TRUE) # 5%
fixest::etable(did_log_r)

etable(list(
  naive_cb_r, twfe_cb_r, did_cb_r,
  naive_one_r, twfe_one_r, did_one_r,
  naive_log_r, twfe_log_r, did_log_r))

etable_output <- fixest::etable(
  list(
    naive_cb_r, twfe_cb_r, did_cb_r, naive_one_r,
    twfe_one_r, did_one_r, naive_log_r, twfe_log_r, did_log_r),
  headers = c(
    "Naïve", "TWFE", "2SDID", "Naïve", "TWFE",
    "2SDID", "Naïve", "TWFE", "2SDID"))

kable_c <- kable(
  etable_output,
  format = "html",
  caption = "",
  escape = FALSE
) %>%
  kable_styling(
    full_width = FALSE,
    position = "center",
    font_size = 14,
    htmltable_class = "classic-table"
  ) %>%
  footnote(
    general = "Significance levels: *** p < 0.001, ** p < 0.01, * p < 0.05, . p < 0.1.",
    general_title = "Note:")

full_html <- paste0(custom_css, as.character(kable_c))
#writeLines(full_html, "kable_c.html")

#Country Clustered Standard Error
naive_cb_r <- feols(cb ~ post_fail_colpus | leaderspell,
  cluster = "gwf_country", data = dt_colpus_r) # 5%
fixest::etable(naive_cb_r)

naive_one_r <- feols(one_cw ~ post_fail_colpus | leaderspell,
  cluster = "gwf_country", data = dt_colpus_r) # 5%
fixest::etable(naive_one_r)

naive_log_r <- feols(cb_log ~ post_fail_colpus | leaderspell,
  cluster = "gwf_country", data = dt_colpus_r) # 5%
fixest::etable(naive_log_r)

twfe_cb_r <- feols(
  cb ~ post_fail_colpus + anticipation +
    t_leader + time_since_last_coup + latent_personalism_lag +
    rgdppc_lag + g_growth_rate_lag + oil_new_lag + per_solider_lag +
    mid_lag + acd_intrastate_lag + exclpop_lag + milethnic_hetero_lag +
    defense_lag | leaderspell + year,
  cluster = "gwf_country",
  data = dt_colpus_r) # 5%
fixest::etable(twfe_cb_r)

twfe_one_r <- feols(
  one_cw ~ post_fail_colpus + anticipation +
    t_leader + time_since_last_coup + latent_personalism_lag +
    rgdppc_lag + g_growth_rate_lag + oil_new_lag + per_solider_lag +
    mid_lag + acd_intrastate_lag + exclpop_lag + milethnic_hetero_lag +
    defense_lag | leaderspell + year,
  cluster = "gwf_country",
  data = dt_colpus_r) # 5%
fixest::etable(twfe_one_r)

twfe_log_r <- feols(
  cb_log ~ post_fail_colpus + anticipation +
    t_leader + time_since_last_coup + latent_personalism_lag +
    rgdppc_lag + g_growth_rate_lag + oil_new_lag + per_solider_lag +
    mid_lag + acd_intrastate_lag + exclpop_lag + milethnic_hetero_lag +
    defense_lag | leaderspell + year,
  cluster = "gwf_country",
  data = dt_colpus_r) # 5%
fixest::etable(twfe_log_r)

did_cb_r <- did2s(
  data = dt_colpus_r,
  yname = "cb",
  treatment = "post_fail_colpus",
  first_stage = ~ t_leader + anticipation + latent_personalism_lag + 
    rgdppc_lag + g_growth_rate_lag + oil_new_lag + per_solider_lag + 
    mid_lag + acd_intrastate_lag +  exclpop_lag + milethnic_hetero_lag + 
    defense_lag | leaderspell + year,
  second_stage = ~ i(post_fail_colpus, ref = 0) +
    t_leader + anticipation + latent_personalism_lag + rgdppc_lag + 
    g_growth_rate_lag +  oil_new_lag + per_solider_lag + mid_lag + 
    acd_intrastate_lag + exclpop_lag + milethnic_hetero_lag + defense_lag,
  cluster_var = "gwf_country",
  bootstrap = FALSE,
  verbose = TRUE) # 5%
fixest::etable(did_cb_r)

did_one_r <- did2s(
  data = dt_colpus_r,
  yname = "one_cw",
  treatment = "post_fail_colpus",
  first_stage = ~ t_leader + anticipation + latent_personalism_lag + 
    rgdppc_lag + g_growth_rate_lag + oil_new_lag + per_solider_lag + 
    mid_lag + acd_intrastate_lag +  exclpop_lag + milethnic_hetero_lag +
    defense_lag | leaderspell + year,
  second_stage = ~ i(post_fail_colpus, ref = 0) +
    t_leader + anticipation + latent_personalism_lag + rgdppc_lag + 
    g_growth_rate_lag +  oil_new_lag + per_solider_lag + mid_lag + 
    acd_intrastate_lag + exclpop_lag + milethnic_hetero_lag + defense_lag,
  cluster_var = "gwf_country",
  bootstrap = FALSE,
  verbose = TRUE) # 5%
fixest::etable(did_one_r)

did_log_r <- did2s(
  data = dt_colpus_r,
  yname = "cb_log",
  treatment = "post_fail_colpus",
  first_stage = ~ t_leader + anticipation + latent_personalism_lag + 
    rgdppc_lag + g_growth_rate_lag + oil_new_lag + per_solider_lag +
    mid_lag + acd_intrastate_lag + exclpop_lag + milethnic_hetero_lag + 
    defense_lag | leaderspell + year,
  second_stage = ~ i(post_fail_colpus, ref = 0) +
    t_leader + anticipation + latent_personalism_lag + rgdppc_lag + 
    g_growth_rate_lag + oil_new_lag + per_solider_lag + mid_lag + 
    acd_intrastate_lag + exclpop_lag + milethnic_hetero_lag + defense_lag,
  cluster_var = "gwf_country",
  bootstrap = FALSE,
  verbose = TRUE) # 5%
fixest::etable(did_log_r)

etable(list(
  naive_cb_r, twfe_cb_r, did_cb_r,
  naive_one_r, twfe_one_r, did_one_r,
  naive_log_r, twfe_log_r, did_log_r))

etable_output <- fixest::etable(
  list(
    naive_cb_r, twfe_cb_r, did_cb_r, naive_one_r,
    twfe_one_r, did_one_r, naive_log_r, twfe_log_r, did_log_r),
  headers = c(
    "Naïve", "TWFE", "2SDID", "Naïve", "TWFE",
    "2SDID", "Naïve", "TWFE", "2SDID"))

kable_cc <- kable(
  etable_output,
  format = "html",
  caption = "",
  escape = FALSE
) %>%
  kable_styling(
    full_width = FALSE,
    position = "center",
    font_size = 14,
    htmltable_class = "classic-table"
  ) %>%
  footnote(
    general = "Significance levels: *** p < 0.001, ** p < 0.01, * p < 0.05, . p < 0.1.",
    general_title = "Note:")

full_html <- paste0(custom_css, as.character(kable_cc))
#writeLines(full_html, "kable_cc.html")

############################# Arellano-Bond GMM ################################

dt_colpus_ab <- panel(dt_colpus, ~ leaderspell + year)
dt_colpus_ab <- dt_colpus_ab %>% rename(id = leaderspell, time = year)
dt_colpus_ab <- pdata.frame(dt_colpus_ab, index = c("id", "time"))

ab_model_cb <- pgmm(
  counterbalancing ~ lag(counterbalancing, 1) + post_fail_colpus + 
    time_since_last_coup + latent_personalism_lag + rgdppc_lag + 
    g_growth_rate_lag + oil_new_lag + mid_lag + acd_intrastate_lag +
    exclpop_lag + milethnic_hetero_lag + defense_lag |
    lag(counterbalancing, 2:99),
  data = dt_colpus_ab,
  effect = "individual",
  model = "twosteps",
  transformation = "d")
summary(ab_model_cb) #Non

ab_model_one <- pgmm(
  one_cw ~ lag(one_cw, 1) + post_fail_colpus + t_leader +
    time_since_last_coup + latent_personalism_lag + rgdppc_lag + 
    g_growth_rate_lag + oil_new_lag + mid_lag + acd_intrastate_lag +
    exclpop_lag + milethnic_hetero_lag + defense_lag |
    lag(one_cw, 2:99),
  data = dt_colpus_ab,
  effect = "individual",
  model = "twosteps",
  transformation = "d")
summary(ab_model_one) #Non

ab_model_log <- pgmm(
  cb_log ~ lag(cb_log, 1) + post_fail_colpus + time_since_last_coup +
    latent_personalism_lag + rgdppc_lag + g_growth_rate_lag +
    oil_new_lag + mid_lag + acd_intrastate_lag +
    exclpop_lag + milethnic_hetero_lag + defense_lag |
    lag(cb_log, 2:99),
  data = dt_colpus_ab,
  effect = "individual",
  model = "twosteps",
  transformation = "d")
summary(ab_model_log) #10%

#Granger Test - Reverse Causality
naive_oneg <- feols(colpus_fail ~ one_cw_lag | leaderspell,
  cluster = "leaderspell", data = dt_colpus)
fixest::etable(naive_oneg) #Non-significant

naive_logg <- feols(colpus_fail ~ cb_log_lag | leaderspell,
  cluster = "leaderspell", data = dt_colpus)
fixest::etable(naive_logg) #Non-significant
```

ALTERNATIVE DEPENDENT VARIABLE

```{r}
#Security Forces Measure
dt_colpus_r <- dt_colpus_r %>%
  mutate(security_force = log(cw_intel + cw_police +
    cw_troops + cw_militia + cw_pg + 1))

sf_naive <- feols(security_force ~ post_fail_colpus | leaderspell,
  cluster = "leaderspell", data = dt_colpus_r) # 5%
fixest::etable(sf_naive)

sf_twfe <- feols(
  security_force ~ post_fail_colpus + t_leader + anticipation +
    latent_personalism_lag + rgdppc_lag + g_growth_rate_lag +
    oil_new_lag + per_solider_lag + mid_lag + acd_intrastate_lag +
    exclpop_lag + milethnic_hetero_lag + defense_lag | leaderspell + year,
  cluster = "leaderspell", data = dt_colpus_r) # 5%
fixest::etable(sf_twfe)

sf_did <- did2s(
  data = dt_colpus_r,
  yname = "security_force",
  treatment = "post_fail_colpus",
  first_stage = ~ t_leader + anticipation + latent_personalism_lag + 
    rgdppc_lag + g_growth_rate_lag + oil_new_lag + per_solider_lag +
    mid_lag + acd_intrastate_lag + exclpop_lag + milethnic_hetero_lag + 
    defense_lag | leaderspell + year,
  second_stage = ~ i(post_fail_colpus, ref = 0) +
    t_leader + anticipation + latent_personalism_lag + rgdppc_lag + 
    g_growth_rate_lag + oil_new_lag + per_solider_lag + mid_lag + 
    acd_intrastate_lag + exclpop_lag + milethnic_hetero_lag + defense_lag,
  cluster_var = "leaderspell",
  bootstrap = FALSE,
  verbose = TRUE)
fixest::etable(sf_did)

etable_output <- fixest::etable(
  list(sf_naive, sf_twfe, sf_did),
  headers = c("Naïve", "TWFE", "2SDID"))

kable_sf <- kable(
  etable_output,
  format = "html",
  caption = "",
  escape = FALSE
) %>%
  kable_styling(
    full_width = FALSE,
    position = "center",
    font_size = 14,
    htmltable_class = "classic-table"
  ) %>%
  footnote(
    general = "Significance levels: *** p < 0.001, ** p < 0.01, * p < 0.05, . p < 0.1.",
    general_title = "Note:")

full_html <- paste0(custom_css, as.character(kable_sf))
#writeLines(full_html, "kable_sf.html")

# Effective Number
dt_colpus_r <- dt_colpus_r %>%
  mutate(effective_number = log(effectivenumber + 1))

ef_naive <- feols(effective_number ~ post_fail_colpus | leaderspell,
  cluster = "leaderspell", data = dt_colpus_r) # 5%
fixest::etable(ef_naive)

ef_twfe <- feols(
  effective_number ~ post_fail_colpus +
    t_leader + anticipation + latent_personalism_lag + rgdppc_lag + g_growth_rate_lag +
    oil_new_lag + per_solider_lag + mid_lag + acd_intrastate_lag +
    exclpop_lag + milethnic_hetero_lag + defense_lag | leaderspell + year,
  cluster = "leaderspell", data = dt_colpus_r)
fixest::etable(ef_twfe)

ef_did <- did2s(
  data = dt_colpus_r,
  yname = "effective_number",
  treatment = "post_fail_colpus",
  first_stage = ~ t_leader + anticipation + latent_personalism_lag + 
    rgdppc_lag + g_growth_rate_lag + oil_new_lag + per_solider_lag +
    mid_lag + acd_intrastate_lag + exclpop_lag + milethnic_hetero_lag +
    defense_lag | leaderspell + year,
  second_stage = ~ i(post_fail_colpus, ref = 0) +
    t_leader + anticipation + latent_personalism_lag + rgdppc_lag +
    g_growth_rate_lag + oil_new_lag + per_solider_lag + mid_lag +
    acd_intrastate_lag + exclpop_lag + milethnic_hetero_lag + defense_lag,
  cluster_var = "leaderspell",
  bootstrap = FALSE,
  verbose = TRUE)
fixest::etable(ef_did)

etable_output <- fixest::etable(
  list(ef_naive, ef_twfe, ef_did),
  headers = c("Naïve", "TWFE", "2SDID"))

kable_ef <- kable(
  etable_output,
  format = "html",
  caption = "",
  escape = FALSE
) %>%
  kable_styling(
    full_width = FALSE,
    position = "center",
    font_size = 14,
    htmltable_class = "classic-table"
  ) %>%
  footnote(
    general = "Significance levels: *** p < 0.001, ** p < 0.01, * p < 0.05, . p < 0.1.",
    general_title = "Note:")

full_html <- paste0(custom_css, as.character(kable_ef))
#writeLines(full_html, "kable_ef.html")
```

APPENDIX: ROBUSTNESS (PT 2011)

```{r}
#Random Forest Imputation for Covariates
pt_vectors <- subset(pt_cw, select = -c(1:12, 45:54))

miceObj <- miceRanger(pt_vectors,
  m = 5,
  returnModels = T,
  verbose = T,
  seed = 777)

pt_vectors_new <- impute(pt_vectors, miceObj, verbose = T)
pt_vectors_final <- pt_vectors_new$imputedData$Dataset_1

#Merge back With Original Variables
pt_dep <- subset(pt_cw, select = c(1:12, 45:54))
pt_cwrf <- as.data.frame(cbind(pt_dep, pt_vectors_final))

#Factor Fixed Effects
pt_cwrf$leaderspell <- as.factor(pt_cwrf$leaderspell)
pt_cwrf$year <- as.factor(pt_cwrf$year)

#Counterbalancing (Sample Average = 1.468)
pt_cwrf <- pt_cwrf %>%
  mutate(cb_count_1 = cbcount + 0.01) %>%
  filter(year %in% c(1946:2010)) %>%
  mutate(mean_cb = mean(cb_count_1, na.rm = T)) %>%
  mutate(counterbalancing = ifelse(cb_count_1 > mean_cb, 1, 0))
mean(pt_cwrf$cb_count_1, na.rm = T) # Sample Average = 1.47

#Log(CB) DV
pt_cwrf <- pt_cwrf %>% mutate(cb_log = log(cbcount + 1))

#Time-To-Treatment
before_after_es_pt <- pt_cwrf %>%
  group_by(leadername, t_leader) %>%
  mutate(ref_yr = case_when(
    coup_fail == 1 ~ t_leader,
    TRUE ~ NA_integer_)) %>%
  ungroup() %>%
  group_by(leadername) %>%
  fill(ref_yr, .direction = "downup") %>%
  mutate(yr_diff = abs(ref_yr - t_leader))

before_pt <- before_after_es_pt %>%
  filter(post_fail == 0) %>%
  mutate(yr_diff = yr_diff - (yr_diff * 2)) %>%
  dplyr::select(c(ccode, year, yr_diff))

before_after_es_pt <- left_join(before_after_es_pt, before_pt,
  by = c("ccode", "year", "leadername"))

before_after_es_pt <- before_after_es_pt %>%
  mutate(duration = if_else(is.na(yr_diff.y), yr_diff.x, yr_diff.y))
before_after_es_pt$duration[is.na(before_after_es_pt$duration)] <- -Inf
before_after_es_pt <- before_after_es_pt %>% arrange(ccode, year)

#Anticipation Effect
before_after_es_pt <- before_after_es_pt %>%
  mutate(
    anticipation = ifelse(duration == -1, 1, 0),
    anticipation2 = ifelse(duration == -2, 1, 0))

#Time Since Last Successful Coup
before_after_es_pt <- before_after_es_pt %>%
  arrange(ccode, year) %>%
  group_by(ccode) %>%
  mutate(
    time_since_last_coup = 0,
    time_since_last_coup = {
      counter <- 0
      sapply(1:n(), function(i) {
        if (coup_success[i] == 1) {
          counter <<- 1
        } else {
          counter <<- counter + 1
        }
        counter
      })
    }
  ) %>%
  ungroup()

# Lag Covariates Within Leader-Spells
dt_pt <- as.data.table(before_after_es_pt)
dt_pt <- dt_pt[order(leaderspell, year)]
dt_pt[, latent_personalism_lag := shift(latent_personalism, n = 1, type = "lag"), by = leaderspell]
dt_pt[, rgdppc_lag := shift(rgdppc, n = 1, type = "lag"), by = leaderspell]
dt_pt[, g_growth_rate_lag := shift(g_growth_rate, n = 1, type = "lag"), by = leaderspell]
dt_pt[, oil_new_lag := shift(oil_new, n = 1, type = "lag"), by = leaderspell]
dt_pt[, per_solider_lag := shift(per_solider, n = 1, type = "lag"), by = leaderspell]
dt_pt[, mid_lag := shift(mid, n = 1, type = "lag"), by = leaderspell]
dt_pt[, acd_intrastate_lag := shift(acd_intrastate, n = 1, type = "lag"), by = leaderspell]
dt_pt[, exclpop_lag := shift(exclpop, n = 1, type = "lag"), by = leaderspell]
dt_pt[, milethnic_hetero_lag := shift(milethnic_hetero, n = 1, type = "lag"), by = leaderspell]
dt_pt[, defense_lag := shift(defense, n = 1, type = "lag"), by = leaderspell]
dt_pt[, post_fail_lag := shift(post_fail, n = 1, type = "lag"), by = leaderspell]

################################ NAIVE OLS #####################################

#All Models Include Lagged IVs With Standard Errors Clustered by Leader Spell
naive_ptcb <- feols(counterbalancing ~ post_fail | leaderspell,
  cluster = "leaderspell", data = dt_pt) #5%
fixest::etable(naive_ptcb)

naive_ptone <- feols(one_cw ~ post_fail | leaderspell,
  cluster = "leaderspell", data = dt_pt) #5%
fixest::etable(naive_ptone)

naive_ptlog <- feols(cb_log ~ post_fail | leaderspell,
  cluster = "leaderspell", data = dt_pt) #5%
fixest::etable(naive_ptlog)

################################## TWFE ########################################

#Full-Specification With Year and Leader-Spell FE (SE Clustered by Leader-Spell)
twfe_ptcb <- feols(
  counterbalancing ~ post_fail + t_leader + anticipation + time_since_last_coup +
    latent_personalism_lag + rgdppc_lag + g_growth_rate_lag +
    oil_new_lag + per_solider_lag + mid_lag + acd_intrastate_lag +
    exclpop_lag + milethnic_hetero_lag + defense_lag | leaderspell + year,
  cluster = "leaderspell", data = dt_pt)
fixest::etable(twfe_ptcb)

twfe_ptone <- feols(
  one_cw ~ post_fail + t_leader + anticipation + time_since_last_coup +
    latent_personalism_lag + rgdppc_lag + g_growth_rate_lag +
    oil_new_lag + per_solider_lag + mid_lag + acd_intrastate_lag +
    exclpop_lag + milethnic_hetero_lag + defense_lag | leaderspell + year,
  cluster = "leaderspell", data = dt_pt)
fixest::etable(twfe_ptone)

twfe_ptlog <- feols(
  log(cbcount + 1) ~ post_fail + t_leader + anticipation + time_since_last_coup +
    latent_personalism_lag + rgdppc_lag + g_growth_rate_lag +
    oil_new_lag + per_solider_lag + mid_lag + acd_intrastate_lag +
    exclpop_lag + milethnic_hetero_lag + defense_lag | leaderspell + year,
  cluster = "leaderspell", data = dt_pt)
fixest::etable(twfe_ptlog)

############################### Two-Stage DiD ###################################

did_ptcb <- did2s(
  data = dt_pt,
  yname = "counterbalancing",
  treatment = "post_fail",
  first_stage = ~ t_leader + anticipation + time_since_last_coup + 
    latent_personalism_lag + rgdppc_lag + g_growth_rate_lag + oil_new_lag + 
    per_solider_lag + mid_lag + acd_intrastate_lag + exclpop_lag + 
    milethnic_hetero_lag + defense_lag | leaderspell + year,
  second_stage = ~ i(post_fail, ref = 0) +
    t_leader + anticipation + time_since_last_coup + latent_personalism_lag + 
    rgdppc_lag +  g_growth_rate_lag + oil_new_lag + per_solider_lag + mid_lag + 
    acd_intrastate_lag + exclpop_lag + milethnic_hetero_lag + defense_lag,
  cluster_var = "leaderspell",
  bootstrap = FALSE,
  verbose = TRUE)
fixest::etable(did_ptcb)

did_ptcb2 <- did2s(
  data = dt_pt,
  yname = "counterbalancing",
  treatment = "post_fail",
  first_stage = ~ t_leader + latent_personalism_lag + rgdppc_lag +
    g_growth_rate_lag +  oil_new_lag + per_solider_lag + mid_lag + 
    acd_intrastate_lag + exclpop_lag + milethnic_hetero_lag + 
    defense_lag | leaderspell + year,
  second_stage = ~ i(post_fail, ref = 0) +
    t_leader + latent_personalism_lag + rgdppc_lag + g_growth_rate_lag +
    oil_new_lag + per_solider_lag + mid_lag + acd_intrastate_lag +
    exclpop_lag + milethnic_hetero_lag + defense_lag,
  cluster_var = "leaderspell",
  bootstrap = FALSE,
  verbose = TRUE)
fixest::etable(did_ptcb2)

did_ptone <- did2s(
  data = dt_pt,
  yname = "one_cw",
  treatment = "post_fail",
  first_stage = ~ t_leader + anticipation + time_since_last_coup + 
    latent_personalism_lag + rgdppc_lag + g_growth_rate_lag + oil_new_lag + 
    per_solider_lag + mid_lag + acd_intrastate_lag + exclpop_lag + 
    milethnic_hetero_lag + defense_lag | leaderspell + year,
  second_stage = ~ i(post_fail, ref = 0) +
    t_leader + anticipation + time_since_last_coup + 
    latent_personalism_lag + rgdppc_lag +  g_growth_rate_lag +
    oil_new_lag + per_solider_lag + mid_lag + acd_intrastate_lag +
    exclpop_lag + milethnic_hetero_lag + defense_lag,
  cluster_var = "leaderspell",
  bootstrap = FALSE,
  verbose = TRUE) #10%
fixest::etable(did_ptone)

did_ptlog <- did2s(
  data = dt_pt,
  yname = "cb_log",
  treatment = "post_fail",
  first_stage = ~ t_leader + anticipation + time_since_last_coup +
    latent_personalism_lag + rgdppc_lag + g_growth_rate_lag + oil_new_lag + 
    per_solider_lag + mid_lag + acd_intrastate_lag + exclpop_lag + 
    milethnic_hetero_lag + defense_lag | leaderspell + year,
  second_stage = ~ i(post_fail, ref = 0) +
    t_leader + time_since_last_coup + anticipation + 
    latent_personalism_lag + rgdppc_lag + g_growth_rate_lag + 
    oil_new_lag + per_solider_lag + mid_lag + acd_intrastate_lag +
    exclpop_lag + milethnic_hetero_lag + defense_lag,
  cluster_var = "leaderspell",
  bootstrap = FALSE,
  verbose = TRUE) #10%
fixest::etable(did_ptlog)
```
